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ABSTRACT 


Spacecraft attitude estimation and pointing aecuracy have always been limited by 
imperfect sensors. The rate gyroseope is one of the most eritieal instruments used in 
spaeeeraft attitude estimation and unfortunately historieal trends show this instrument 
degrades signifieantly with time. Degraded rate gyroseopes have impaeted the missions 
for several NASA and ESA spaeeeraft, ineluding the Hubble Teleseope. A possible 
solution to this problem is using a mathematieally modeled dynamie gyroseope in lieu of 
a real one. In this thesis, data from sueh a gyro is presented and integrated into a 
spaeeeraft attitude estimation algorithm. 

The impediment to spaeeeraft attitude estimation presented by imperfect sensors 
has been overeome by developing more aeeurate sensors and using Kalman filters to 
reduee the effeet of noisy measurements. Kalman filters for spaeeeraft attitude 
estimation have historieally been based on an Euler angle or quaternion formulation. 
Though Euler angles and quaternions are arguably the most eornmon methods with whieh 
to deseribe the attitude of a spaeeeraft, other methods of deseribing attitudes do exist - 
ineluding the Gibbs and Rodriguez parameters. A Kalman filter based upon the Gibbs 
parameter is presented and analyzed in this thesis. 
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I. INTRODUCTION 


A. SPACECRAFT ATTITUDE ESTIMATION 

The pointing accuracy requirements of modem Department of Defense satellites 
have increased steadily over the past few decades. As the pointing accuracy requirements 
for a satellite increase, so does the level of accuracy in the estimation of its attitude. 
Spacecraft attitude estimation is an extremely complex and non-linear process. Inputs 
from several different sensors - such as star trackers and rate gyroscopes, among others - 
are required to determine the attitude of a spacecraft. Because no sensor is error-free - 
even at the time of manufacture - and all man-made equipment degrades with age, the 
problem of accurate attitude estimation throughout the mission life of a satellite is 
extremely important. A method for overcoming these impediments via a Kalman 
filtering process is presented and analyzed here. 

Kalman filtering has been used in spacecraft attitude estimation for quite some 
time. The earliest published application was in 1970 by Farrell and several others have 
followed since. Lefferts, Shuster, and Markley published a thorough review of the topic 
in 1982 and since then Markley, Crassidis, and several others have kept Kalman filtering 
an active topic of research in the space industry [Markley]. Several different attitude 
representations have been used in Kalman filtering with varying degrees of success. 
Palermo successfully implemented a Kalman filter using an Euler angle representation of 
the attitude for a simulated bifocal relay mirror spacecraft [Palermo] and his dynamic 
model is used as a starting point for this work. 

B, SATELLITE ENVIRONMENT 

The environment in which a satellite operates complicates the problem of 
spacecraft attitude estimation and control. A satellite in orbit is subject to non-constant 
external torques at all times - some secular (varying linearly with time), some periodic, 
and some random. The forces that induce these torques are gravity gradient, magnetic, 
solar pressure, and atmospheric drag, though the effects of atmospheric drag on satellites 
outside of LEO is considered negligible. Though these torques can be predicted with 
high accuracy (with the exception of atmospheric drag), any errors in prediction couple 
with errors in sensor accuracy to increase the error in attitude estimation. Because large 
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errors in attitude ean be attributed to the effeets of these external torques, a brief 
overview of gravity gradient, magnetie, and solar torque and the errors that ean result 
from eaeh is presented here. For more in depth information the reader is referred to texts 
by Sidi, Bong Wie, or Hughes [Sidi, Bong Wie, Hughes]. In the rest of this seetion we 
survey a number of disturbanees aeting on the spaeeeraft. 

1. Gravity Gradient Torque 

A spaeeeraft is not a point mass and may not be treated as sueh. It is a rigid body 
(in most oases) with a mass distribution about a oenter of mass. Vioe treating it as a point 
mass an inertia dyadio is used. Gravity gradient torques are imparted to a spacecraft 
because gravity acts on each element of the spacecraft. Gravity acting on a mass m 
located at a distance r from the spacecraft center of mass will induce a torque about the 
center of mass. The effects of gravity will act on each portion of the satellite in 
accordance with Newton’s Laws of Gravitation. Ignoring the effects of the moon and 
other third bodies due to their small effects, these torques may be written as: 

( 1 - 1 ) 

where ju is the gravitational parameter for the Earth, is the distance from the satellite 
to the center of the Earth, and is the inertia dyadic for the satellite. 

Some satellites with low pointing requirements actually use gravity gradient 
stabilization. A quick examination of equation 1-1 reveals that knowledge of the position 
of the satellite and the inertia dyadic are critical for predicting the gravitational torque. 
System identification algorithms are being actively researched and developed to facilitate 
more precise knowledge of spacecraft inertia dyadics. 

In addition to its own mass distribution, the fact that the gravitational field 
produced by the Earth is an aspherical potential further complicates the problem. Eike a 
satellite, the Earth does not have a uniform mass distribution - hence its aspherical 
gravity field. Zonal, tesseral, and sectoral harmonics within the field are both non-linear 
and extremely complex to model mathematically. Inaccuracies in modeling the 
gravitational field of the Earth can lead to attitude pointing errors as well. 
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2, Magnetic Torque 

The Earth has a rotating magnetic field B through which any orbiting satellite 


must travel. The satellite itself has an intrinsic magnetic moment m = 




. Interaction 


between the magnetic moment of the spacecraft and the magnetic field of the Earth 
generates a torque on the spacecraft calculated by 


Ng =mxB 


( 1 - 2 ) 


The magnitude of 5 will decrease as the satellite altitude increases. Eor dealing 
with torques induced by the interaction of the spacecraft and Earth magnetic fields, 
precise knowledge of the spacecraft magnetic moment m is critical. As with the moment 
of inertia, system identification algorithms can be used to mitigate errors due to 
inaccurate estimates of fh . 

3. Solar Torque 

Maxwell’s equations imply that electromagnetic waves have momentum, which 
may be transferred to objects with which it comes in contact. Since light is an 
electromagnetic wave, it exerts pressure. Though this pressure is miniscule in an Earth 
environment it is not miniscule for a satellite in orbit about the Earth. 

Eor electromagnetic radiation, basic physics shows us ihai Electromagnetic Force 
= Work + Energy Density which, when expressed in a more mathematical manner 
becomes 

\ {Tm)da = \ FdV + —\ ^dV (1-3) 

js Jv dt^^c 

where T is the Maxwell stress tensor, F is the force density, and S is the Poynting 

vector {S = -^ExB). 

471 

The amount of pressure exerted on an object by an electromagnetic wave is highly 
dependent upon the type of surface being illuminated. In some complex models, the type 
of reflection - specular or diffuse - as well as the reflectivity of the surface plays a part in 
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the calculation. In a simple model the following equation is often used to determine the 
maximum solar torque imparted to a spacecraft 


N. 


SP 


F 

= — 4 (1 + ^) cos(/)(c - c ) 
c 


(1-4) 


where is the spacecraft area, q is the reflectance factor, i is the angle of 
incidence of the incoming light, -c^ is the distance between the spacecraft center of 
gravity and the center of solar pressure, and is the solar power density (which varies 
with time). Solar cycles, changes in -c^ due to fuel expenditures, and the changing 

area of a spacecraft tracking a point on the Earth all contribute to the difficulty of 
estimating and compensating for disturbances due to solar pressure [SMAD]. 

4, Atmospheric Drag 

As previously mentioned the effects of atmospheric drag on satellites outside of 
LEO are considered negligible. Eor LEO satellites, however, atmospheric drag is the 
most difficult external torque to predict. Because the dynamics of the outer reaches of 
the atmosphere are not fully understood it cannot be modeled accurately. The effects of 
drag on LEO spacecraft are directly proportional to the area of the spacecraft. 
Atmospheric drag, though pertinent to LEO applications, was not included in any of the 
models used in developing this thesis but is mentioned here for completeness. 

C. SENSORS 

Accurate measurement of data from external sources is required for a spacecraft 
attitude control system to estimate its attitude. Though several sensors exist that perform 
this function - including sun sensors, horizon sensors, and earth sensors among others - 
the star tracker is the most accurate and pertinent to the focus of this work. 

1. Rate Gyroscopes 

A cursory perusal of either the Euler Equations of Motion or the kinematic 
equations for the quaternion or Gibbs parameters reveals their dependence upon angular 
rate data. Any error in angular rate measurement or calculation will result in an error in 
spacecraft attitude estimation. Angular rate information is critical to the accurate 
estimation of spacecraft attitude - regardless of the estimation method. Rate gyroscopes 
provide this information to the spacecraft attitude control system. 
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There are several different kinds of rate gyroseopes available for use in the spaee 
environment today. Meohanieal rate gyroseopes and rate integrating gyroseopes have 
been in use on orbit for quite several years. Both of these types of gyroseopes are 
dependent upon meohanieal moving parts whioh degrade over time. Laser gyroseopes, 
quartz rate sensors, and hemispherioal resonator gyroseopes are now available today with 
muoh higher aoouraoies and no relianoe on moving parts [Sidi]. 

2, Dynamic Gyroscopes 

Rate gyroseopes are man-made devioes and therefore have finite lifetimes - 
espeoially in the harsh environs of spaee. Beoause all of the external seoular and periodio 
torques on a spaoeoraft may be modeled with some modieum of aeeuraey and the torques 
applied to the spaeeeraft via the momentum exehange deviees of its attitude eontrol 
system are known with a fair amount of precision, it is possible to determine the angular 
rate of the spacecraft via mathematical modeling and the use of external measurements. 
Such an algorithm is called a dynamic gyroscope. The use of a dynamic gyroscope upon 
failure of a simulated mechanical gyroscope is shown and analyzed in this work. 

3, Star Trackers 

Navigators from ancient times used the stars as navigational aids. Satellites 
navigating in space do the same via a device called a star tracker. Because stars may be 
considered inertially fixed bodies for all intents and purposes and because they are 
extremely small as seen from our solar system, they are ideal objects to use as an attitude 
reference. Star trackers provide the most accurate attitude data to the spacecraft attitude 
control system by several orders of magnitude over any other type of sensor. While older 
star trackers were capable of tracking only one star at a time, the new generation of star 
trackers can feed attitude quaternions to its host satellite and track multiple stars 
simultaneously. A new type of star sensor is currently being developed by Dr. Junkins et 
al which will use two simultaneous images of star fields from different parts of the sky to 
determine spacecraft attitude [Junkins]. 

D, KALMAN FILTERS 

As previously mentioned, a satellite attitude control system receives input from 
imperfect external sensors to estimate the actual attitude of the satellite. The Kalman 
filter is an estimation algorithm frequently employed to do this. In essence, the Kalman 
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filter is a set of mathematieal equations that provides a reeursive solution of the least- 
squares estimation method. It supports eomputational estimations of past, present, and 
future states and ean do so even when the nature of the system under consideration is not 
precisely known. 

Given the extremely non-linear and varying nature of the space environment, this 
makes it an ideal tool with which to perform spacecraft attitude estimation. The 
effectiveness of the filtering algorithm is highly dependent upon how the state vector is 
defined, the type of filter employed, and several other factors. In this research we 
developed a Kalman Filter based on specific parameterization of the spacecraft attitude 
and attitude error. 

E, STAR GAPS AND RATE GYROSCOPE UPSETS 

Satellites in orbit currently experience periods in which their star trackers are 
unable to sense stars - called star gaps. Rate gyroscopes in some older satellites are 
sometimes sending data over 1000 times the actual reading - such an event is called a 
rate gyroscope upset. Both of these events wreak havoc upon attitude control algorithms 
- the Kalman filter in particular. Both of these phenomena are simulated and their effects 
mitigated via adaptive covariance and dynamic gyroscope integration. 

F. THESIS ORGANIZATION 

In this thesis we introduce the dynamics of spacecraft attitude in Chapter Two. 
Chapter Three begins with a cursory overview of the Kalman filtering process followed 
by a brief section on attitude error representation. It closes with a derivation of the Gibbs 
parameter based Kalman filter. Chapter Four covers the problems caused by star gaps 
and methods of managing such problems. Chapter Five introduces the dynamic 
gyroscope and methods of dealing with rate gyroscope upsets. 
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II. SPACECRAFT ATTITUDE 


A. ATTITUDE PARAMATRIZATON 

In this section we review spacecraft attitude estimation methods. There are 
several methods for doing this, but only a short description of most of them will be 
included here. Emphasis will be placed on the quaternion because understanding the 
quaternion is critical to comprehending the development of the Kalman filters used in 
spacecraft attitude estimation. 

Regardless of the method, describing the attitude of a spacecraft basically entails 
describing the orientation of one coordinate system with respect to another. For purposes 
of satellites, a coordinate system based on the principal axes of the spacecraft - called the 
body coordinate system - and Earth Centered Inertial (ECI) systems are used. 

1. Direction Cosine Matrix 

The direction cosine matrix is the simplest manner in which to describe the 
attitude of a spacecraft. Given two coordinate systems each consisting of three 
orthogonal unit vectors there exists a three by three matrix C that relates the two 
coordinate systems. 

When dealing with different coordinate systems it becomes necessary to develop 
a notation for annotating in which coordinate system a vector is expressed. In this paper 
a superscript to the left of the vector will indicate the coordinate system in which a vector 
is expressed. For example, 'v is a vector in the inertial reference frame while *v is a 
vector in the body reference frame. 


A vector x = 


expressed in body coordinates may be expressed in inertial 


coordinates by multiplying it by a direction cosine matrix (DCM) as shown below 

‘C^'’x = ‘x ( 2 - 1 ) 


Note the two superscripts on the DCM. For a direction cosine matrix, the 
superscript on the left indicates the coordinate system the DCM transforms a vector to 
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while the right superseript indieates the eoordinate system from whieh the veetor is being 
transformed. Conversely, a veetor in the inertial frame may be expressed in body frame 
eoordinates by 

'^C f = '’x (2-2) 

As one might expeet 'C* and are related as the transpose of eaeh other. 

2, Euler Angles 

Euler angles are a series of rotations about the axes of one eoordinate system to 
align it to another eoordinate system. These rotations are sometimes referred to as roll cp 
(rotation about the x axis), piteh 0 (rotation about the y axis), and yaw y/ (rotation about 
the z axis). 

Sinee in the Euler angle formulation the order matters, there are twelve possible 
eombinations of Euler angles one may use to go from one eoordinate system to another. 
Multiplieation of the individual rotation matriees for eaeh Euler angle will result in a 
direetion eosine matrix. Given the two eoordinate systems a and b the direetion eosine 
matrix for a roll-piteh-yaw Euler sequenee would be 

10 0 eos^ 0 -sin^ eos^^ sin^^ 0 

“C* = 0 coscp sin^ 0 10 -sin^ cos^ 0 (2-3) 

0 -sin^ eos^3 sin^ 0 cos^ 0 0 1 

3, Quaternions 

A mathematieal strueture ealled the quaternion is a eonvenient method with whieh 
to deseribe the orientation of a eoordinate system. 

a. Quaternion Definition 

Quaternions are a form of hyper-eomplex numbers invented by William 
Hamilton in the nineteenth eentury and are today used extensively both in roboties and 
spaeeeraft attitude eontrol. Quaternions are represented as a four element set eonsisting 
of three veetor eomponents and one sealar eomponent 
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There is no set convention for which element of the quaternion is a scalar, 
but in this paper will always be the scalar element of the quaternion. 

b. Attitude Quaternions 

When used for describing an attitude, quaternions are constrained to the 
surface of a four dimensional hyper-sphere defined by 

+ql+ql+ql=\ (2-5) 

To determine the physical meaning of quaternions constrained on this 
hyper-sphere, consider two coordinate systems as shown below, S and S’. The vector 
portion of the quaternion represents an axis about which coordinate system S must be 
rotated to align it with the 5” coordinate system - called the eigenaxis. 

s 



Figure 1. Coordinate Systems S and S’ 






The scalar portion of the quaternion is a measure of the magnitude of the 


rotation 6*^. Defining the eigenaxis as a unit vector e = 


be defined as follows: 


= sm 

= sin 
^3 = sin 
^4 = cos 


U j 

U j 
(oA 

U J' 
ro 

U j 


, the quaternion elements may 


(2-6) 


Since quaternions are vectors of unit magnitude, they may not be added 
together through the standard definition of addition. Quaternion algebra is a rich and 
interesting topic, but one that will not be dealt with here in its entirety. Two uses of 
quaternion algebra will be addressed because of their usefulness in attitude estimation - 
the quaternion conjugate and quaternion multiplication. 

c. The Quaternion Conjugate 

The quaternion conjugate of quaternion q is given by 

q* =-q^-q^-q^+q^ ( 2 - 7 ) 


From a physical standpoint, the conjugate of a quaternion q represents a 
rotation of the same magnitude about a vector in the opposite direction. It is also worthy 
to note that the inverse of a quaternion is identical to its complex conjugate. The 
quaternion conjugate is useful when using the quaternion as a rotation operator, which 
will be discussed later. 

d. Quaternion Multiplication 

As noted previously, quaternions have their own unique algebra. The 
symbol 0 is used to denote quaternion multiplication. It can be shown that the 
quaternion product is defined by 
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p®q = 


dA 

di 

-dl 

dl 

Pi 


dA 

dl 

dl 

Pi 

dl 

-dl 

dA 

dl 

Pi 

_-dl 

-dl 

-dl 

dA_ 



( 2 - 8 ) 


which may also be written as 


pQ9q = 


p^q + q.p-pxq 

PSa-P^ 


(2-9) 


The quaternion product is a useful tool for time propagation of the 
quaternion - as will be shown later. Recall that a quaternion represents a rotation from a 
reference frame to a given attitude. It naturally follows that the quaternion product 
defined in equation 2-9 would also represent a rotation - whieh it does. Defining A{q) as 
a DCM representing the same rotation as quaternion q, it may be shown that 

A{p®q) = A{p)A{q) (2-10) 

e. Quaternion Propagation 

Given an angular veloeity vector o) in body coordinates and an initial 
orientation expressed by the quaternion q^^ the orientation at any time t may be expressed 
as the solution of the following differential equation 


dq 

dt 


= q = 



( 2 - 11 ) 


with the initial eondition ^(0) = q^,. 

f. Advantages of Quaternions 

Quaternion representations of spacecraft attitude hold several advantages 
over the Euler angle and direction cosine matrix representations. Perhaps the most 
obvious advantage comes from the size - a quaternion is a four element structure, 
whereas a direction cosine matrix has nine elements. A quaternion describes the same 
attitude with half of the number of elements, saving both memory and proeessing power. 
The simplieity of the time derivative of the quaternion makes it an ideal method with 
which to do attitude propagation. Also of great importance is the faet that it has no 
singularities. 
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4, Gibbs Vector 

Gibbs vectors are closely related to quaternions, though they exist in □ ^ viceD . 
The Gibbs vector is defined as 


^ 1 
g = — 

^4 


^2 

% 


= e tan 


UJ 



( 2 - 12 ) 


The Gibbs vector is obviously singular whenever goes to zero, which occurs 
for rotations of 180 degrees. Consequently, the Gibbs parameter is an extremely 
effective tool for describing rotations in the interval (-180°,+180°). The relationship 
between the Gibbs vector and the quaternion is shown in the gnomonic projection in 
figure 2; 



Figure 2. Gibbs vector as a Gnomonic Projection 

B. ATTITUDE SENSORS 

Measurement of external data is required in order to produce an attitude estimate. 
There are a multitude of sensors available to measure attitude data for a spacecraft, 
including horizon sensors, sun sensors (coarse and fine), and star trackers. None of these 
provides a perfect measurement, hence the need for attitude error estimation. A brief 
overview of rate gyroscopes and star trackers along with their respective sources of error 
is presented here, as both are simulated in models used in this thesis. 
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1, Rate Gyroscopes 

Rate gyroscopes provide angular rate data to the spacecraft attitude control 
system. There are several different types of rate gyroscopes available today, including 
mechanical, fiber optic, laser, and hemispherical resonator gyroscopes. Regardless of the 
type of gyroscope, all perform the same function. Each sends the spacecraft attitude 
control system the angular rates of the satellite in body coordinates. It is beyond the 
scope of this work to describe each in detail, but one thing that all rate gyroscopes have in 
common is that they are inherently noisy. The errors sent to the attitude control system 
by noisy rate gyroscopes may be categorized into four distinct types of error; bias error, 
random walk, scale factor error, and orthagonality errors. 

a. Bias 

Bias error is sometimes referred to as a constant drift. Bias errors for 
current generation gyroscopes range from 0.00 T / hour to T / hour. The measured 
angular rate for a rate gyroscope may be expressed as 

= (2-13) 

With b(t) representing the bias at time t. As will be shown later in this 
work, it is possible to estimate gyroscope bias via a Kalman filtering process. 

b. Rate Walk 

Rate walk is sometimes described as bias drift. Random walk propagates 
from the white noise of the sensor and is responsible for non-deterministic behavior. Its 

units are normally "/yjhour . 

c. Scale Factor 

The scale factor error is the linear deviation of the measured rate from the 
true rate - normally given as a percentage or in parts per million. Asymmetry and non¬ 
linearity have been observed in scale factors for rate gyroscopes. Scale factors are caused 
by imperfections in manufacturing and the degradation of the gyroscope with the passage 
of time [Hewitson et al]. Though a source of error in any gyroscope, scale factors are not 
used in any of the models in this work, but are presented here for completeness. 
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d. Orthogonality Errors 

Rate gyroscopes are mounted in clusters relative to one another to provide 
data to a spacecraft attitude control system. Any misalignment of the sensors with 
respect to one another will result in error - known as orthogonality error. Launch 
vibrations and thermal deformation over time could cause such errors to occur. As with 
scale factor errors, no orthogonality errors are modeled in this work - they are presented 
for theoretical completeness. 

2, Star Trackers 

As their name implies, star trackers generate attitude data based upon the relative 
position of stars as seen from the spacecraft itself. Each star sensor has a star catalog 
with which to compare the image it sees. Due to the large number of recorded stars, the 
exact composition of the catalog will be determined by the orbit in which the host 
satellite will be placed - there is no such thing as a standard star catalog. 

a. Star Characteristics 

In order to use stars for navigational purposes it is necessary to 
differentiate between stars. There are two characteristics of stars which make this 
possible - magnitude and spectra. 

Magnitude refers to the brightness of a star. There are two types of 
magnitude - absolute and apparent. The apparent magnitude is how bright an object 
seems when viewed from the Earth. In contrast, the absolute magnitude is the apparent 
magnitude of an object placed at ten parsecs away. Star trackers use the apparent 
magnitude for their calculations. 

The spectrum of a star refers to the type of radiation it is emitting. Stellar 
spectra are divided into seven categories. There are O, B, A, E, G, K, and M class stars, 
with O being the hottest and M being the coolest star. 

b. Star Sensing 

Current generation star trackers use Charge Coupled Devices (CCDs) to 
detect stars, though future generation star trackers will use CMOS technology [Junkins]. 
A stray light shield prevents light from outside the star tracker bore sight from reaching 
the CCD. Once the CCD has registered a star field image, the star sensor processor will 
match that image with its library and send an attitude quaternion to the spacecraft attitude 
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control system. This is an extremely simplified explanation of what actually occurs, but 
the aetual meehanies of this proeess are outside the seope of this work. 

c. Star Tracking 

Onee one or more stars have been acquired, the spaeecraft will move 
around the Earth in its orbit and possibly change its attitude intentionally. The software 
must track the image of the star through this motion and be ready to aequire a new star 
when the image exits the Field of View (FOV) of the star traeker [Sidi]. 

d. Star Gaps 

There are times when a star exits the field of view (FOV) of a star tracker 
and no new star within the new FOV has been aequired. During sueh time periods no 
data is being sent to the satellite attitude eontrol system - such periods are ealled star 
gaps. Star gaps may also be caused by the satellite pointing in such a way that the moon 
or sun prevents the star traeker from sensing any stars due to their high intensity. Star 
gaps are detrimental to attitude determination when they extend for more than a few 
seeonds and can wreak havoe on attitude estimation algorithms, as will be shown in a 
later seetion. 

e. Star Tracker Error 

Although easily the most aecurate sensor available, star trackers are not 
error free. Fike rate gyroscopes, star trackers also have a scale factor, which will vary 
with age. The CCD assemblies are temperature sensitive structures - a primary source of 
error. Misalignments of star trackers with respect to the spacecraft body as a result of 
launch vibration or thermal deformation over time are also a souree of star traeker error. 
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III. KALMAN FILTERING FOR ATTITUDE ESTIMATION 


As shown in the preceding section, satellite attitude control systems estimate the 
attitude of the satellite based upon measurements from imperfect sensors. As time passes 
the sensors degrade and the data they send to the attitude control system degrades as well. 
The Kalman filter provides a solution to this problem. In this section the Kalman filter is 
introduced in its basic form. Attitude error representations are introduced and then a state 
vector is defined for attitude estimation purposes. A discrete Kalman filter for attitude 
estimation based upon the Gibbs parameter is then developed and the results of its 
implementation are presented. 

A, THE KALMAN FILTER 

As alluded to previously, a Kalman filter is a recursive mathematical algorithm 
that allows one to estimate the state of the system based upon its previous states, its 
known dynamics, and knowledge of the accuracies of the sensors involved. Once 
initialized, there are three steps to the filtering process: prediction, measurement, and 
correction. A diagram of these steps is shown below. There are two main types of 
Kalman filters - regular and extended. A brief overview of each is presented here. 


Measurement 




t 

Initialization 

Figure 3. Kalman Filtering Process Diagram 


Correction 
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1, The Kalman Filter 

The Kalman filter may be used to prediet the future state of a proeess governed by 
the following equation 

^k^,=Fk^k+Gu, + w, (3-1) 

where k represents the time step and x is a veetor of state variables x e 91”. is 
an nxn matrix known as the state transition matrix which relates the current state to the 
next one. The matrix G is a nxl matrix that relates the control input at time step k- - 

to the state x. The final term is the process or plant noise. This noise accounts for 
random inputs and inaccuracies of the state space model that cause the plant to perform in 
a non-deterministic manner. The noise is assumed to be white noise with a normal 
probability distribution 

p{w)UN{Q,Q) (3-2) 

where is a matrix representing the plant noise w defined by 

Q,=cov(w,) (3-3) 

Although the state vector represents all the variables being estimated, not all the 
estimated variables are required to be measured. The measurement is related to the state 
by the following equation 

^k=HkX,+v, (3-4) 

whereis an mxn matrix often referred to as the measurement matrix or the 

measurement sensitivity matrix in some cases. The second termv^ represents the 

measurement noise - or sensor inaccuracies. As with the plant noise, it is assumed to be 
white noise with a normal distribution given by 

p{v,)U (3-5) 
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where is the measurement error eovarianee matrix defined by 

R,=co\(v,) (3-6) 

It is important to note that though they share a similar probability distribution, the 
plant and measurement noise are eompletely independent of one another. If a system 
under eonsideration ean be deseribed by equations 3-1 and 3-3, then a diserete Kalman 
filtering algorithm may be used to estimate its future states. 

a. Prediction 

The first phase of a filtering algorithm is the predietion phase. In this 
phase at time step k a predietion is made for the next time step based upon the 
eurrent state using the system model. This predietion is given by 

(3-7) 

the minus in the superseript of x denotes that the estimate has been made 
prior to taking a measurement. 

The error eovarianee is also predieted during the predietion stage. The 
eovarianee matrix P is a u x n matrix representing the estimate error eovarianee. The 
eovarianee matrix is predieted via 

P,.,=F,PX + Q, (3-8) 

b. Measurement 

Onee a predietion for the next state has been made, an aetual measurement 
is taken. Note onee again that it is not neeessary to measure every state being 
estimated. A quantity known as the residual is then ealeulated by 

Zk=Zk-Hi^x^ (3-9) 

A quantity ealled the Kalman gain is then defined aeeording to 
Kk =PkHl{HkPkHl +7?,)-' (3-10) 
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where 


Pk =E\ix,-x,){x,-xJ] 

and 

(3-11) 

^k =E[x,\z^_^] 

K=E[x^\z^'\ 

c. Correction 

The Kalman gain in eonjunetion with the measurement residual are then 
used to find the estimate via 


(3-12) 

The Kalman gain is also used to update the eovarianee matrix by 
P:=iI-K,H,)P;{I-K,H,)y+K,R,Kl (3-13) 

Onee this has been eompleted, the proeess starts over again - henee the 

term reeursive. 

d. Filter Initialization 

In order to begin the Kalman filtering proeess an initial estimate of both 
the state and the eovarianee matrix are required to begin the reeursion. The filter ean be 
initialized in several different ways, but the most common way to do so is via a simple 
guess. The closer to the actual state the initial guess is, the faster the estimate generated 
by the filter will reach the actual value. On occasion the filter is initialized using the first 
few measurements, but since it is easier to simply make an educated guess this method is 
not often used. 

2, The Extended Kalman Filter 

Not all processes that are desired to be estimated can be modeled in the linear 
form required for the Kalman filter. For these processes, the Extended Kalman Filter 
(EKE) can be used. The EKE is simply a regular Kalman filter that linearizes about the 
estimate and covariance. 

For the EKE, the state vector x is governed by a non-linear function/dependent 
upon the current state, a control input, and plant noise 
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(3-14) 


with the measurement z being expressed as a non-linear funetion of the state, as 

Zk=h{x,,v^) (3-15) 

Both noise terms and have the same properties as with the linear Kalman 

filter. Sinee it is impossible to know the noise at any given time and both the plant and 
measurement noise are characterized as white noise, both the state and measurement 
vectors are approximated by setting the noise equal to zero 

Zk=KXk,0) (3-16) 

Linearizing about the zero mean non-linear functions in equation 3-16 will yield a 
close approximation to the actual value of the function itself - if the noise is truly zero- 
mean in nature, its expected value must be zero by definition. We linearize the system by 
taking the Jacobians of f and h with respect to the state vector , and the following 
matrices are obtained 


T7 

(3-17) 

dx, 

(3-18) 


With the above matrices defined the EKF algorithm proceeds in exactly the same 
manner as the regular Kalman filter. 

a. Prediction 

The prediction for the EKF is done via the non-linear function/described 
previously. Eor the predicted state and covariance the equations are 

K+x=fi^k^u„0) (3-19) 

Pk.^=FkPkFl^Qk (3-20) 
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b. Measurement 

Once a prediction has been made a measurement is taken and as 
before, a residual z^. is calculated, this time via the non-linear function h 

^k=^k-hiXk,0) (3-21) 

The Kalman gain is calculated as before in equation 3-9, with the 
exception that when using the EKF H is the Jacobian of the non-linear function h vice the 
sensitivity matrix //used in the regular Kalman filter. 

c. Correction and Initialization 

As with the regular Kalman filter, the EKF corrects the estimate and the 
covariance matrix based upon the measurement and the Kalman gain. The equations 
used to do this are identical to the regular Kalman filter. Also like the regular Kalman 
filter, an EKE must be initialized with a starting value - like the correction step of the 
filtering process, the EKE method for initialization is identical to that used for the regular 
Kalman filter. 

B, ATTITUDE ERROR REPRESENTATIONS 

When using Kalman filters for spacecraft attitude estimation, the state vectors 
most frequently consist of bias and attitude errors. Prior to developing a Kalman filter for 
attitude estimation, it is appropriate to cover how the attitude error may be represented. 
The two attitude error representations relevant to the discussion here are the quaternion 
error and the Gibbs parameter error representations. Emphasis will be placed upon the 
latter of the two. 

1, Quaternion Error 

If two reference frames are slightly offset from one another, the error quaternion 
Sq represents the rotation that will align one frame with the other. Eor example, if the 
estimated attitude of a spacecraft is given by and the measured attitude of the 
spacecraft by a perfect sensor was ^ then Sq would represent the rotation , the 

estimated attitude, to^, the actual attitude. Eor all rotations it can be shown that 

q = dq®q^^f (3-22) 
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The ease of proeessing the above equation makes the error quaternion a useful 
representation for spacecraft attitude error. However, from a Kalman filtering standpoint, 
the error quaternion presents some problems. Recall that when defining a residual, z^. 

the estimated value is subtracted from the measured value. For illustrative purposes, 
define p as the measured error quaternion and q as the predicted error quaternion based 
upon the current state. Subtracting to obtain the residual yields 

a 1 r^ii Ta-^i 

Pi A A-A 

z=p-q= 

A A A-A 

_aJ LaJ La-a 

While this seems straightforward enough, recall that attitude quaternions reside on 
the surface of a hyper-sphere defined previously and shown again here 

ql +ql+ql +ql = l (2-5) 

Though p and q both satisfy the requirements of equation 2-5 their difference will 
not do so, and therefore is not an attitude quaternion. This becomes a problem when 
using quaternions in a Kalman filtering algorithm, so converting the quaternion to 
different representations becomes necessary - as will be shown later in this chapter. 

2, Gibbs Error 

Although the Gibbs parameter is closely related to the quaternion - as one can 
easily see from its definition in chapter two - it does not have the normalization 
constraint. The Gibbs parameter for a particular quaternion resides in a plane tangent to 
the surface of the hyper-sphere on which quaternions reside. A direct mapping between 
the Gibbs parameter plane and the quaternion hyper-sphere exists and is given by 



Sq{ag) = {A + a^) ^ 



(3-24) 


It can further be shown that the direction cosine matrix associated with a 
particular 5q{a) may be approximated by 

-[ax]-^(aV 3 , 3 -aa^) (3-25) 
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This direct relationship between the Gibbs parameter and the quaternion will be 
used in the next section to formulate a Kalman filter for spacecraft attitude estimation 
[Markley]. 

C. FORMULATION OF A GIBBS PARAMETER BASED KALMAN FILTER 

With the Kalman filter and attitude error representation introduced, the 
formulation of the Kalman filter may be developed. The Kalman filter developed in this 
section estimates the spacecraft attitude error and the gyroscopic bias. It uses a six 
element state vector of the form 


a{t) 

bit) 


(3-26) 


where a(t) represents the attitude error in Gibbs parameters and b(t) is the 
gyroscopic bias error. 

1, Filter Overview 

The Kalman fdter developed here is based on one developed by F. Landis 
Markley [Markley]. It treats the actual attitude as the quaternion product of the estimated 
quaternion error and the previous estimate, as shown here 

qit) = dqiait)) 0 q^^^ (3-27) 


where q^^^- is the attitude estimate from the previous time step and Sqiait)) 
represents the rotation from the last time step to the attitude at the current time step. 

The steps for the filter are the same as for any other - prediction, measurement 
update, and correction. What is unique about this filter is that the propagation of the state 
ait) 


vector X = 


bit) 


is relatively trivial in nature, as Markley himself states. Prior to 


developing the filter itself, it is of interest to show why this is the case. 

In the prediction step of this filter, the attitude quaternion from the previous time 
step ^^y^is considered the optimal estimate for the attitude of the spacecraft, i.e. the 

attitude error estimate a (t) at the beginning of every time step is zero by definition. 
Once a measurement is taken, is no longer the optimal estimate for the spacecraft 
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attitude. During the eorreetion step, the eorreeted error estimate a^{t) is used to 
propagate the old estimate, so that it beeomes the new optimal estimate for the 
spaeeeraft attitude via 

q:.f=5q{a\t))®q;^f (3-28) 

Onee the new optimal estimate has been generated, d{t) is immediately reset to 
zero. When the next step in the recursion begins, the new then becomes the optimal 
estimate 

q^ref = Sqia (0) ® q^ref = ^^(0) ® q^ref = q^ref (3-29) 

Because the filter is designed in this manner, the actual propagation of the state 
vector is trivial because its first three elements are always reset to zero. The value of 
immediately after the correction step of the filtering process is the actual attitude 

estimate - the state vector merely tracks the error. 

2, System Dynamics - Developing the State Transition Matrix 

The first step in formulating any Kalman filter is mathematically modeling how 
the system should behave and subsequently constructing a state transition matrix based 
upon this mathematical model. For a satellite, we have its current state quaternion and its 
angular rate o). From equation 2-12, it is known that the time derivative of the quaternion 
is based upon the angular rate of the body it is describing. 

Placing this in terms of the filter formulation 

9.„=i ®?„/ (3-30) 

where is the angular rate of the satellite at the reference attitude, q^^^. 

The filter itself is not estimating the quaternion components of the attitude - it 
estimates the error in terms of Gibbs parameters. The equation for the propagation of the 
Gibbs parameter error is given by 

1 r 1 

= (Ax 3 + - ) - - (^+«../) X «g = /(«. 0 
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(3-31) 



as shown by Markley [Markley]. 

We now propagate the estimate from time nT to time nT + T . Therefore 

^ref - = eonstant (3-32) 

and = 0 sinee constant. Therefore, from equation 3-27 we obtain 
q{nT + T) = dq{a{nT + T))®q{nT) (3-33) 

Propagating in a discrete manner, a^{nT + T) may be written as 
a^{nT + T) = a^{nT) +d^T (3-34) 

When this approach is taken, equation 3-31 may be rewritten as 

1 r 1 

ag=il3^3+-agag)(i)--(i)xag (3-35) 


Linearizing about yields the following equation 


d 


g 


1 

CO — coxa„ 
2 ^ 


(3-36) 


Recalling from equation 2-16 that the angular velocity is actually a function of the 
measured velocity and the gyroscopic bias and substituting that into equation 3-36 
(assuming a zero noise component for the time being) gives 


d = CO —b — — (co —b)xa 

g meas ^ ^ ^ ' g 

— CO -b-—Q. a 

meas ^ S 


(3-37) 


Introducing the term which is defined as the skew of the measured angular 
velocity 


Q = X = 

m L rneas 




-co^ 










(3-38) 



We propagate the attitude from time nT to time nT+T equation 3-34 may be 
rewritten as 


a^{nT + T) = a{nT) + {cD^^JnT)-b{nT)-^QJnT)a{nT))T 
= ih.,-^ClJnT)T)a(nT) + co^^^inT)T-binT)T 


(3-39) 


Assuming a constant bias over short periods of time T, the state equation 3-39 
may be used to formulate the state transition matrix for the Kalman filter 


4 = 


/ Q T -I T 

-*3x3 2 “ 3x3-^ 


0 


3x3 


‘3x3 


(3-40) 


which will propagate the state as 

x{nT + T) = AjX{nT) + u{nT) + w{nT) (3-41) 

This is in the standard state space form so we are able to apply standard Kalman 
filtering algorithms with the control input modeled as 

u = coT (3-42) 

3. System Sensors - Developing the Measurement Sensitivity Matrix 

With the state transition matrix developed from the system dynamics, it becomes 
necessary to determine how the filter will take an external measurement and incorporate 
it into the filtering algorithm via the measurement sensitivity matrix. The term which 
drives the ‘correction’ equation is the error between the actual measurement and the 
predicted measurement. In particular in our problem we have the actual orientation 'v 
from the star trackers library in the inertial coordinate frame. At the same time, given an 
estimate of the spacecraft orientation q and the observed star *v from the star tracker we 
can compute and estimate 

‘v = A(q)^v (3-43) 

In particular let 

q = dq{a)®q^{nT) (3-44) 
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so the error ‘v - ‘v ean be written as 

'v-'v= 'v - A{dq{a))A{q^ {nT)Yv (3-45) 

Substituting equation 3-35 into the above equation produees a relationship 
between the Gibbs error a, q ^^^, *v and 'v 

'v =[43 -[«x]-^(a4,3 (3-46) 

In the above equation q^^j. is the estimate from the previous time step. 

Condueting a first order expansion with respeet to the Gibbs error about yields the 
following equation 

- Bh - - 

hCv) = hCv) + — [ax] 4 = hCv) + H^a (3-47) 

5 V 


From this, is shown to be 

Bh 

4-^ ['vx] (3-48) 

With the above relations developed, the measurement sensitivity matrix ean be 
formed by relating them to the state veetor x 


// = [4 03,3] (3-49) 

4. The Filter Equations 

With the state transition and measurement sensitivity equations eonstrueted, it is 
now possible to summarize and implement the Kalman filter. In this version of the filter, 
the predietion and eorreetion steps are eombined for expedieney. 

The Kalman filter is built on the following state spaee model, developed above 


x{n +1) = 2 “ x{n) + + w = Ax Bu + w 

0 0 l43_ 

L *^3x3 *^3x3 J 

~ Bh 

y(n)= [*vx] 03 ^ 3 ] x{n) + v = Hx + v 

B V 


(3-50) 
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where x{n) = 


and u{n) = . 


a 

b 


The following equations are used to estimate the state and the covariance 


= 




‘3x3 ^ “m 


0 


3x3 


3x3-* 


‘3x3 


Xk^K{z^-{{ vx] 03 + M 


(3-51) 


= AjX,+K{z^-Hx^) + u 


Pk.r = AA4 + Q-Kk{HP,H^+R)Kj (3-52) 


Note that the above two equations differ slightly from the Kalman filtering 
equations derived earlier. For implementation purposes, the prediction and correction 
steps were combined. It is also important to note here that the measurement is the 
measured error, which is given by 


z 


k 




(3-53) 


The Kalman gain K used in equation 3-60 is calculated by 

K,=A,P,H\HP,H^ +Ry (3-54) 

A diagram of the filtering process is shown in Figure 4. The reference quaternion 
which is updated after the filter estimates the attitude error gives the optimal estimate of 
the attitude of the spacecraft after it is updated with the information from the filter. Once 
the reference quaternion becomes the optimal estimate, the attitude error estimate is reset 
to zero. 
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Figure 4. Gibbs Parameter Kalman Filter Flow Diagram 


5. Testing the Filter 

Onee developed, the filter was implemented in MATLAB and tested using both 
‘elean’ and eorrupted quaternion measurements taken from a simulated satellite that 
drives itself to nadir pointing after an initial offset. The gyroseopes in this simulation had 
no bias. Clean quaternions and gyroseope readings were used on this initial run to see if 
the Kalman filter would funetion eorrectly. The results are shown in figures five and six. 
Though some quaternion error is present, it is on the order of 10and is attributed to 
numerieal error only. The estimated bias is of the same order and is also attributed to 
numerieal errors and the fact that a non-zero value for the bias was used to initialize the 
filter. Because of the low level of error, it was concluded that the filter functions 
correctly. 
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Quaternion Errors 
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Figure 5. Quaternion Errors For Gibbs Parameter Kalman Filter With Zero Bias and No 

Measurement Noise 
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IV. STAR GAP ERROR MITIGATION 


As a satellite travels in its orbit about the Earth, there will be times when none of 
its star trackers are able to feed an attitude quaternion into the attitude control system. 
There are several possible reasons for these star gaps, including solar and lunar entry into 
the field of view (FOV) of the star tracker, processing time for the tracker to find and 
lock on a new star(s) once the previous one(s) exit the FOV, or due to a reorientation 
maneuver by the satellite - among others. These gaps have marked effects upon attitude 
estimation algorithms. In this section the effects of star gaps are examined - both for the 
Gibbs parameter based Kalman filter developed in Chapter Three and an Euler angle 
based Kalman filter developed by Palermo using the same data [Palermo]. Methods to 
mitigate the effects of these star gaps are presented along with simulated results of their 
implementation. 

A, EFFECTS OF STAR GAPS 

When a star gap occurs, the effect it has on the attitude estimation algorithm is 
highly dependent upon the formation of that algorithm. Regardless of the type of 
algorithm in use, a longer star gap equates to a larger error as will be shown in the results 
later in this section. 

1, Euler Angle Based Kalman Filter 

Prior to examining how star gaps affect the Euler angle based Kalman filter it is 
useful to briefly examine its formulation. For a full derivation of the equations, the 
reader is referred to [Palermo]. When a star gap occurs in the Euler angle based Kalman 
filter, the attitude estimate diverges for all four elements of the quaternion. Upon 
reacquisition of a star, there is a brief spike in the attitude error then the measurement 
updates drive the error back down - eventually. The results of a 200 second star gap on 
an attitude estimator based upon Euler angles are shown in the following figure. The 
estimator is tracking noisy measurements for both attitude and angular rate. 
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Figure 7. Euler Angle Based Kalman Filter with 200 Seeond Star Gap 

2, Gibbs Parameter Based Kalman Filter 

The Gibbs Parameter based Kalman Filter behaves in a very different manner. 
When a star gap oeeurs, the estimator is able to continue the track with excellent accuracy 
because it is still receiving angular rate data. The measured angular rate data allows the 
estimator to propagate the reference quaternion with little loss of accuracy due to the 

star gap. The graph was made using the exact same data for a 200 second star gap as was 
used for the Euler angle based Kalman filter. Eigure 8 shows the quaternion errors for 
the same track. 
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Q1 Track vs. Actual 



Q3 Track vs. Actual 



Q2 Track vs. Actual 



Q4 Track vs. Actual 



Figure 8. Estimated Quaternion Elements vs. Actual for Gibbs Parameter Based Kalman 

Filter (Estimate Shown in Blue) with 200 Second Star Gap 


Quaternion Errors - One Hz Sample Rate 



Eigure 9. Quaternion Error for Gibbs Parameter Based Kalman Eilter with 200 Second Star 

Gap 
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B, ADAPTIVE COVARIANCE 

One technique for mitigating the effects of a star gap is adapting the different 
covariance matrices intrinsic to the Kalman filtering process. The three matrices dealing 
with covariance are the filter covariance matrix P, the plant covariance matrix Q, and the 
measurement covariance matrix R. The P matrix is predicted and corrected based upon 
measured values and the Kalman gain - adapting it based upon situational parameters 
would significantly change the way in which the filter propagates the estimate. The plant 
and measurement covariance matrices, however, remain constant throughout the filtering 
process. The performance of any Kalman filter may be modified by changing their value. 
By doing so adaptively based upon external conditions, it is possible to mitigate some of 
the attitude error caused by intermittent star gaps. Adaptation of the plant and 
measurement covariance are examined individually and then examined in conjunction 
with one another. 

C. THE PLANT COVARIANCE MATRIX 

The plant covariance Q may be thought of as a measure of how well the model 
emulates the actual dynamics of the quantities being estimated. The more precise the 
model, the lower the values contained in the Q matrix. A star gap might be 
conceptualized as a system modeled in an extremely poor manner, i.e. a system with an 
extremely large plant noise. Using this conceptualization, a trigger was developed that 
increased the plant covariance when a star gap occurred. This covariance trigger was 
implemented on both the Euler and Gibbs parameter based Kalman Filters. 

1. Adapting Plant Noise for the Euler Angle Based Kalman Filter 

The Euler angle based Kalman Filter used for testing the adaptive covariance 
concept was developed by Palermo using SIMULINK [Palermo]. This format was 
maintained, but modified to adapt the covariance matrices. The following figure shows 
the trigger that was developed for this purpose. The figure shown is the measurement 
covariance trigger - the plant noise covariance trigger is identical. An initial matrix 
serves as the covariance trigger until a certain point - user determined - at which it is 
increased by a user determined factor. After a certain time period has elapsed, the 
covariance matrix is reset to its initial value. A digital clock was used to do this, as can 
be seen in the diagram. 
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Figure 10. Measurement Covarianee Trigger 

Onee the eovarianee trigger was implemented, it was tested using several different 
gains for the 200 seeond star gap under eonsideration. Note that the trigger was timed to 
inerease the eovarianee just before the star gap and deerease it immediately just after the 
star gap - assuming that the measurements preeeding the star gap were stored in memory 
and reeyeled through with the new value of Q onee the star gap was deteeted. Both the 
magnitude and the direetions of the rotations varied for eaeh gain used. The gains for the 
plant noise eovarianee matrix were inereased logarithmieally for the duration of the star 
gap. The following figures shows the mean quaternion errors plotted against the gains 
applied to the plant noise matrix Q for the duration of the star gap. 
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Figure 11. Mean Quaternion Error for Plant Noise Injection During Star Gaps 

2, Adapting Plant Noise for the Gibbs Parameter Based Kalman Filter 

The Kalman filter based on Gibbs Parameters was built in MATLAB without 
using SIMULINK, making any adaptation of the Q matrix easy to accomplish. Unlike 
the Euler angle based Kalman filter, adapting the plant noise with the Gibbs Parameter 
based Kalman filter produced no noticeable effect. 

D. THE MEASUREMENT COVARIANCE MATRIX 

The measurement covariance matrix R provides the filter with the kn own 
accuracies of the sensors being used to take the measurements z. Conceptually, it follows 
that when a sensor is producing no measurement its function should come into question. 
This was done for both types of Kalman filters discussed in this paper. 

1, Perturbation of Measurement Noise Covariance Matrix for an Euler 
Angle Based Kalman Filter 

Adapting the measurement covariance matrix accordingly using the trigger shown 
in figure eight produced the following results when applied to the Euler angle based 
Kalman filter. 
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Figure 12. Mean Quaternion Error for Measurement Matrix Adaptation During Star Gap 

2, Perturbation of Measurement Noise Covariance Matrix for a Gibbs 
Parameter Based Kalman Filter 

As with the plant noise injection, measurement noise injection for the Gibbs 
parameter based Kalman filter produced no difference with regards to estimation error 
encountered during a star gap. 

E. SIMULTANEOUS ADAPTATION OF COVARIANCE MATRICES 

If adjusting the sensitivity of the filter to plant and measurement noise 
individually affects the average quaternion error for a Kalman filter, it follows that 
adapting both simultaneously would do so as well. This was accomplished via use of two 
triggers operating in parallel, as shown previously in figure eight. The results are shown 
in the following figure - a significant improvement over adapting either matrix 
individually. 
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A\«rage Quaternion Error with Plant and Measurement Noise Injection 


O 



10 10 
Plant and Measurement Matrices Gains 


Figure 13. Average Quaternion Error with Simultaneous Plant and Measurement Covarianee 

Matrix Adaptation 

Adaptation of the plant noise and measurement noise matriees is an option for 
mitigating attitude estimation error in some eases, as has been shown in the preeeding 
seetions. While it does not eliminate the error indueed due to star gaps, it does have the 
potential to reduee the amount of error eneountered. 
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V. DYNAMIC GYROSCOPE AND RATE GYRO UPSETS 


Mechanical rate gyroscopes have historically been a weak point in satellite 
design. Their high failure rate coupled with their noisy output and importance in attitude 
estimation makes them a liability from a reliability standpoint. Improved gyroscopes 
have been developed - such as the laser gyroscope - but many systems currently in orbit 
are still relying on mechanical gyroscopes which are both noisy and prone to failure. 

This chapter introduces the dynamic gyroscope - a software alternative to the rate 
gyroscope, shows how rate gyroscope upsets can affect attitude estimation algorithms, 
and then integrates the dynamic gyroscope with attitude estimation algorithms to mitigate 
the effects of rate gyroscope upsets. 

A, THE RATE GYROSCOPE UPSET 

The rate gyroscope upset is a phenomenon that has been occurring on certain 
satellites in orbit. At certain random points in time the rate gyroscope will send a 
measurement to the attitude control system that is 100 to 1000 times the actual reading. 
This wreaks havoc on the control system, as can be seen in the following two figures. 
Both figures were produced with the same data from a simulated satellite - with a 
simulated rate gyroscope upset. The simulated rate gyroscope upset was of two seconds 
duration with a measurement 1000 times the correct reading being sent to the attitude 
estimation system. As can be seen by looking at the figures, the Gibbs Parameter based 
Kalman Filter developed earlier performs extremely well as compared to the Euler angle 
based filter. After a brief jump in estimator error, the attitude error for the Gibbs 
Parameter Kalman filter returns to the order oflO "^ while the Euler Angle Kalman filter is 
unable to recover any semblance of accuracy after an upset occurs. 
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Quaternion Error - Euler Angle Model with Two Second Rate Gyroscope Upset 
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Figure 14. Quaternion Error with 2 second Rate Gyroscope Upset - Euler Angle Based 

Kalman Eilter 


Quaternion Error - Gibbs Parameter Model with Two Second Rate Gyroscope Upset 



Eigure 15. Quaternion Error with 2 second Rate Gyroscope Upset - Gibbs Parameter Based 

Kalman Eilter 
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Rate gyroscope upsets of longer and longer lengths were tested on the Gibbs 
Parameter Kalman filter, as shown in figure 16 and up to a certain point the attitude 
estimator was able to recover its pre-upset estimation accuracy. After a certain time 
period the estimator does lose its ability to recover from an upset, as shown in figure 17 
when the rate gyroscope upset was 150 seconds in length. However, its superiority over 
the Euler Angle Kalman filter remains apparent. 


Quaternion Error - Gibbs Parameter Model with 95 Second Rate Gyroscope Upset 



Figure 16. Gibbs Parameter Based Kalman Filter with 95 second Rate Gyroscope Upset 
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Quaternion Error - Gibbs Parameter Model with 100 Second Rate Gyroscope Upset 


- q1 Error 

- q2 Error 


•1 _ - q3 Error 

- q4 Error 



_2 L_[_I_1_^_1_[_I_1_I__ 

0 50 100 150 200 250 300 350 400 450 500 

Time 


Figure 17. Gibbs Parameter Based Kalman Filter with a 150 Seeond Rate Gyroscope Upset 

A bit of analysis revealed that the covariance matrix for the Kalman Filter was 
positive semi-definite for all in which it was able to recover from the rate gyroscope 
upset. If the covariance matrix became non-positive semi-definite during the rate 
gyroscope upset, the fdter would not recover. 

B. THE DYNAMIC GYROSCOPE 

It is possible to mathematically model and predict all of the external torques on a 
satellite with some degree of accuracy - some more than others. Recall from basic 
physics that the torque A of a satellite during a time period may be expressed as 

«—=['■ A* (5-1) 

•'4 

Equation 5-1 may further be broken down because the torques on a satellite may 
be expressed as the sum of its external torques - which may be mathematically modeled 
- and its internal torques from any momentum exchange devices. Equation 5-1 may 
therefore be written as 
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+ (5-2) 

•’h 

It follows that 

iV„„ = kdt = -H,^ (5-3) 

Assuming that any change in the moment of inertia between and 4 will be 
negligible and reealling that 

(5-4) 

It follows that 

^4., ~^h =^K+i-4) 

(5-5) 

4+1 = 4 +^ ‘(^4. -^0 

Inserting equation 5-3 into equation 5-5 gives a way to predict the angular 
momentum for 4 ^j given only known or mathematically modeled data 

4+1 “ 4 ext ^MEd) ( 5 " 6 ) 

Using equation 5-6 or some variant thereof, the angular momentum of a satellite 
may be obtained analytically. An algorithm that calculates the angular rate in this 
manner is referred to as a dynamie or a pseudo-gyroscope. Such a gyroscope was 
developed by Aerospace Corporation and a variant of this was implemented by Palermo 
[Palermo]. The following figure shows a comparison of a simulated real gyroscope and 
the Palermo dynamic gyroscope. The dynamie gyroscope produees a reading very close 
to that produced by the real gyroscope. While rate error is not absent, the amount of error 
- shown in figure 19 - is such that the resulting attitude error will be manageable. 
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Figure 18. Comparison of Real and Dynamic Gyroscope Performance 


Dynamic Gyroscope - Gyroscope Difference 



Figure 19. Difference Between Real and Dynamic Gyroscope Readings 
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Though convenient, dynamic gyroscopes are not free of problems. Any 
inaccuracy in knowledge of the moment of inertia of the spacecraft will result in angular 
rate errors, as can be seen from equation 5-6. Any expenditure of fuel, movement of fuel 
due to any slew maneuver or satellite motion, or any movement of the solar arrays - 
among other things - will change the moment of inertia of a satellite, thereby eausing 
some error. System identifieation algorithms ean be used to help mitigate this problem, 
but sueh algorithms are not perfect - error will always remain. 

C. INTEGRATING THE DYNAMIC GYROSCOPE WITH AN ATTITUDE 

ESTIMATOR 

By using a dynamic gyroscope in eonjunction with a real meehanical gyroseope it 
is possible to mitigate the effect of a rate-gyroscope upset. As seen in the above figures, 
the dynamie gyroscope is able to produee data quite close to that of a real gyroseope. 
Given the fact that the Gibbs Parameter based Kalman filter is extremely robust with 
regards to recovering from rate gyroscope upsets and the faet that any mechanism for 
switching from regular gyroscope readings to measurements from the dynamie gyroscope 
will be imperfect, the Gibbs Parameter based Kalman filter is an ideal estimator with 
whieh to integrate the dynamic gyroscope. The following figure shows the quaternion 
error for a Gibbs parameter based Kalman filter with a 95 seeond rate gyroscope upset of 
1000 times the actual reading. The dynamie gyroseope in this seenario took over 
immediately so no upset readings entered the estimator. 
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Figure 20. Comparison of Estimator Performance with and without Dynamic Gyroscope 

input during a 95 second Rate Gyroscope Upset 

As can be seen, when a dynamic gyroscope is in use and no data from the real 
gyroscope during its upset period reaches the estimator, the amount of error encountered 
is only slightly increased. As mentioned earlier, this would be an ideal scenario. It 
would be impossible to design a perfect algorithm where the dynamic gyroscope data 
would replace the real gyroscope measurements immediately. Some upset data would 
always have a chance of getting into the estimator. The robustness of the Gibbs 
parameter based Kalman filter to short rate gyroscope upsets was shown previously. The 
following figure shows a scenario in which the dynamic gyroscope takes over two 
measurements after the commencement of the rate gyroscope upset. Notice the minor 
perturbations in the error - the estimator rapidly converges back to an estimate with error 
on the order of 10^. 
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Quaternion Errors - Rate Gyroscope Upset with Imperfect Timing of DG Takeover 



Figure 21. Quaternion Error with Dynamie Gyroscope Takeover Occurring Two Seconds 

after Rate Gyroscope Upset Commencement 

By creating a switching algorithm with a high reliability, rate gyroscope upsets 
may be experienced by an attitude estimator using a Gibbs parameter based Kalman filter 
with little loss of estimator accuracy, as can be seen in the figures above. 
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VI. CONCLUSIONS 


The Gibbs parameter formulation of the Kalman filter provides an aeeurate 
estimate of spaeeeraft attitude. It is far more robust an estimator when eonfronted with 
either a star gap or rate gyroseope upset than the Euler angle estimator. 

A, SUMMARY 

A Kalman filter based upon the Gibbs parameterization of spaeeeraft attitude is 
developed and analyzed. A eomparison between this formulation and a Kalman filter 
developed by Palermo is eondueted under different operating seenarios. The dynamie 
gyroscope is shown to be a viable substitute for mechanical gyroscopes during periods of 
rate gyroscope upsets. The Gibbs parameter based Kalman filter used in conjunction 
with the dynamic gyroscope is shown to be an excellent method to mitigate the effects of 
rate gyroscope upsets and produces far superior results than the Euler angle formulation 
of the Kalman filter when faced with a rate gyroscope upset. 

B, RECOMMENDATIONS 

It is strongly recommended that future work be done on attitude estimation in two 
specific areas at a minimum: Kalman Eiltering and Unscented Eiltering. 

1. Kalman Filtering 

The Kalman filter developed in this work was limited to six states for purposes of 
the research. It is possible - and somewhat desirable in some cases - to have a much 
larger state vector estimating other factors impacting spacecraft attitude estimation such 
as scale factors of gyroscopes and star trackers, alignment errors in the sensor to body 
transformation matrices, misalignment of gyroscopes, and errors in the ECI to body 
coordinate transformations - among others. Grey, Kolve, Herman, and Westerlund have 
developed such a filter that is used in on-orbit calibrations [Gray]. The Aerospace 
Corporation has developed a similar filter. It is recommended that such filters and the 
modeling to support their implementation be developed and studied - implementing the 
Gibbs parameter formulation for attitude error estimation into these models. 

2, Unscented Filtering 

As mentioned in the description of the EKE, it is based upon a linearization of 
non-linear functions. This linearization is the fundamental flaw of the EKE. Julier et al 
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have developed an alternative to the EKF based on the premise that it should be easier to 
approximate a Gaussian distribution than to approximate an arbitrary non-linear function. 
Markley and Crassidis have adapted this filtering method to spacecraft attitude estimation 
and it is recommended that research into this method as an alternative to Kalman filtering 
be conducted [Crassidis]. 
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